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Thermo-ylscoolaatlc  Stressea  In  a  Sphere 
with  an  Ablating  Cavity*. 

T.  G.  Rogers  and  E.  H.  Lee 
Brov/n  University 

Abstract 

Thermo-vlscoelastlc  stresses  In  a  sphere  with  a  con¬ 
centric  spherical  cavity  are  analysed,  for  uniform  initial 
temperature,  and  cavity  ablation  at  an  elevated  constant  surface 
temperature.  The  material  Is  considered  to  be  thermo -rheological 
ly  simple,  with  properties  prescribed  by  an  isothermal  relaxation 
modulus  function  and  a  relaxation  time  scale  factor  given  as  a 
function  of  the  temperature.  The  theory  is  developed  in  the  form 
of  integral  equations  in  time  for  the  radial  stress  gradient . 
These  lend  themselves  to  nximerlcal  solution,  and  an  example  la 
presented  of  a  sphere  of  polymethyl -methacrylate  for  which  the 
relaxation  modulus  function  and  temperature -log (time)  shift  func¬ 
tion  are  available.  Features  of  the  solution  are  discussed,  and 
compared  with  a  corresponding  elastic  case  and  the  previously 
published  solution  for  an  ablating  solid  sphere. 


This  work  was  sponsored  by  the  Bureau  of  Weapons,  Department 
of  the  Navy,  under  Contract  NOrd  1859^  with  Brovm  University. 
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’  Introduction  ^ 

Tills  paper  Is  concerned  with  thermal  stresses  in  a 
linear  viscoelastic  material  in  which  increase  in  temperature 
produces  thermal  expansion  with  coefficient  o,  and  also  contracts 
the  time  scale  of  the  relaxation  modulus  (and  creep  compliance) 
by  a  factor  9(T),  where  T  Is  the  temperature  measured  above  a 
base  temperature  T^.  The  response  to  stress  at  different 
temperatures  of  many  materials  can  be  expressed  by  relations  of 
this  type  [l]*,  termed  thermo-rheologically  simple  behavior  by 
Schv/arzl  and  Staverman  [2],  and  known  in  the  chemical  literature 
as  the  Williams -Landel -Perry  law.  On  the  common  log(tlme)  plot, 
this  corresponds  to  change  in  temperature  causing  a  shift  of  the 
relaxation  curve  along  the  log (time)  axis  without  change  of  shape 
In  order  to  deal  with  varying  (time-dependent)  tempera¬ 
ture,  with  the  resulting  contraction  of  the  time  scale  thus 
Itself  a  function  of  time,  it  is  convenient  [3]  to  define  a 
reduced  time  K  which  Incorporates  the  varying  scale  factor  t?(T), 
so  that  in  tenns  of  ^  the  viscoelastic  stress -strain  relation  at 
varying  temperature  has  the  same  form  as  the  isothermal  law  at 
the  base  temperature  T  ; 

K  ^e.  .(x,  ^') 

3ij(x,?)  =  J  2G(^-^«)  ,  (1) 

“  o 

where  x  denotes  the  triplet  of  rectangular  Cartesian  coordinates 

•» - 

Numbers  in  square  brackets  refer  to  the  bibllogi'aphy  at  the 
end  of  the  paper. 
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s. .  and  e . .  are  the  stress  and  strain  deviators 
X  d  3  IJ  IJ 

respectively,  and  G(t)  is  the  relaxation  modulus  in  shear  at  the 
base  temperature  The  reduced  time  ^(x,t),  where  t  Is  the 

real  time,  is  determined  in  terms  of  the  varying  scale  factor 
9[T(x,t)]  by  the  relation 

t 

5(x,t)  =  J  9[T(x,t*)]dt’.  (2) 

o 

During  the  course  of  the  analysis  v^e  shall  find  it 
convenient  to  consider  a  dependent  variable,  for  example  the 
stress  devlator,  sometimes  to  be  a  function  of  the  reduced  time 
5:s^j(x,5),  and  sometimes  expressed  in  terms  of  the  real  time, 
and  we  shall  use  the  same  function  notation:  s^j(x,t).  This 
signifies  replacing  ^  in  s^^Cx,?)  by  ^(x,t).  The  functional  form 
is  thus  different,  but  no  confusion  need  arise  since  the  independ¬ 
ent  variables  will  either  be  stated  or  be  clear  from  the  context . 

^s .  ^ 

Space  derivatives,  e.g.  will  always  signify  differentiation 

at  constant  real  time  t,  since  they  v;ill  arise  only  from  the 
equilibrium  and  displacement-strain  relations,  which  involve 
space  gradients  at  fixed  real  time.  Since  a  transformation  from 
t  to  ^  maj'  be  involved  in,  for  example,  the  interpretation  of  an 
Integral  over  reduced  time  of  a  space  derivative,  a  distinguishing 
function  notation  would  become  cumbersome  in  the  present  problem. 

The  integral  operator  form  of  the  viscoelastic  stress - 
strain  I'elation  (l)  is  utilized  because  of  its  generality 


Only  the  relationship  (l)  for  shear  type  deformation  has  been 
stated  since  we  shall  assume  the  response  in  dilatation' to  be 
elastic , 
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compared  v/ith  the  lovr-order  differential  operator  laws  commonly 
considered  in  viscoelastic  stress  analysis  solutions.  This 
generality  permits  accurate  representation  of  material  response 
over  the  equivalent  of  a  wide  frequence  band,  and  is  needed  in 
thermal -stress  problems  because  of  the  sensitivity  of  the  scale 
factor  9  to  temperature  change.  Moreover,  the  integral  operator 
law  can  lead  to  representation  of  the  stress  distribution  sought 
as  the  solution  of  an  Integral  equation  with  its  kernel  dependent 
on  the  relaxation  modulus  G(t) .  Numerical  solution  by  finite 
sum  methods  may  prove  convenient  [4],  and  permits  direct  substitu 
tion  of  the  experimentally  measured  material  characteristics  into 
the  analysis. 

In  this  paper  we  treat  the  problem  of  a  sphere  with  a 
concentric  spherical  cavity.  The  temperature  is  initially 
uniform  throughout,  when  that  at  the  cavity  surface  increases 
discontlnuously  to  a  value  v/hlch  is  thereafter  maintained 
constant,  and  the  cavity  ablates  with  steadily  increasing  radius 
until  it  reaches  the  outer  radius,  and  the  viscoelastic  material 
has  then  disappeared.  Such  a  problem  falls  within  the  general 
framework  of  the  method  devised  by  Mukl  and  Sternberg  [5]»  but 
the  moving  boundary  rules  out  application  of  the  Laplace  trans- 
fonn,  and  Integral  operator  expressions  must  be  maintained  to 
yield  integral  equation  formulations  which  can  cope  with  the 
moving  boundary  [4].  This  aspect  of  the  analysis  was  developed 
in  [6],  but  it  was  applied  there  only  to  the  problem  of  a  solid 
sphere  in  which  the  coupling  between  the  elemental  spherical 
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shells  happens  effectively  to  disappear  because  of  symmetry 
requirements  at  the  center.  In  the  present  problem  the  spherical 
shell  solutions  are  linked  by  the  stress  boundary  conditions  at 
the  surfaces,  and  this  leads  to  coupling  in  the  space  variables 
of  the  Integral  equations  in  time,  and  so  to  a  much  more 
challenging  computational  task.  This  coupling,  incidentally, 
rules  out  the  transform  approach  even  for  fixed  boundaries, 
since  the  resulting  integral  equation  is,  in  general, -no  longer 
of  convolution  type,  as  pointed  out  by  Sternberg  and  Gurtin  [7]- 

The  spherically  symmetric  stress  field. 

A  body  in  the  form  of  a  spherical  shell  is  considered, 
loaded  xiniformly  over  each  surface  and  subjected  to  a  radially 
symmetrical  temperatvire  distribution  T(r,t),  where  r  is  the 
radius.  A  brief  development  of  the  theory  given  in  [6]  is 
presented  here  to  provide  a  frcunework  for  discussion  of  the  new 
solution.  Symmetry  determines  the  principal  stresses  to  be  the  . 
radial  stress  d  and  the  circumferential  stress  repeated.  The 
sum  of  the  principal  stresses  is  therefore: 

<j  =  d^+2d^.  (3) 

The  single  non-trivlally  satisfied  equilibrium  equation  in 

Symmeti’y  determines  a  single  independent  stress  deviator  compo¬ 
nent,  and  a  corresponding  strain  devlator  component,  which  for 
convenience  can  be  replaced  by  principal  stress  and  strain 
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differences.  V/ith  u(r,t)  as  the  radial  displacement,  the 
equivalent  of  (l)  then  becomes 

Og(r.K}  -  =  f  20(K-K’)  ir  If  -  gjJdS'  (5) 

o 

v/hlch,  using  (4),  can  be  >frltten  in  the  form 

.  (6) 

Note  that,  as  mentioned  In  the  Introduction,  the  r  derivatives 
operate  on  functions  of  r  and  t,  but  transformation  to  r  and  K 
Is  Implied  for  the  Integration,  dg  can  be  eliminated  from  (3) 
and  (4)  to  give: 

Ij  -  A  .  (7) 

Elastic  dilatatlonal  response  is  assumed,  together  with  thermal 
expansion  with  constant  coefficient  a,  giving  the  relation 

0  =  3k[|a  .1.  2m  _  3„0],  (8) 

where  k  is  the  bulk  modulus  and  0  the  Increase  in  temperature 
over  the  uniform  initial  temperature.  Integration  of  (7)  at 
constant  real  time  and  substitution  for  0  from  (8)  yields 

r^cij^(r,t)-[a(t)  ]^d^[a(t),t]  =  3k[r^u(r, t) -[a( t)  ]^u[a(t) , t  ] 

r  p 

-  3a  J  .P  6(p,t)dp]  ,  (9) 

a(t) 

in  vfhich  conditions  at  a  boundary  radius  a(t)  give  the  lower 
limits  of  the  integrals.  Substitution  for  u(r,t)  from  (9)  into 
(6)  gives: 
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P  -St 

3/= -4  J  o(^-c‘)  5^  +  ^( 


3  ,cJj,(a,t»)  u(a,t‘' 


^  C)  II  ■'■  a%(a)}]d?«  (10) 

T*  • 


r  a 


The  integral  Is  a  convolution  integral  (see,  for  example, 
Tricomi  [8]),  and  this  provides  a  simplification  of  (lO)  through 
the  associative  property  of  convolutions .  We  define  the  function 
R(0  by  the  Integral  equation 


R(^)  +  ^  J  G(^-^«)  dC'  =  2G(C) 


(11) 


which,  it  v/ill  be  noted,  contains  the  same  operator  on  R  as  does 
3<j 

(10)  on  3^.  The  associative  property  of  convolutions  then 
enables  (lO)  to  be  replaced  by: 


=  i  R(^-?')|^  f- 


3  .d„(a,t «) 


u(a,t • . 


-  ^  {  J  If  <3p  +  a-^©(a)}].  (12) 

r  a  ^ 

We  consider  an  ablating  hollovf  sphere,  and  take  a(t)  to  be  the 
varying  cavity  radius.  The  surface  boundary  conditions  are  then: 


r  =  a(t),  c5^(a(t),t]  =  f3^(t);  r  =  b(t),  cJ^[b(t),t]  =  f2(t)  (13) 

v/here  bCt)  is  the  radius  of  the  external  surface,  and 
-f-j^(t),  -f2(t)  are  prescribed  surface  pressures.  Since  u(a, t) 
is. not  known  a  priori,  the  Integrand  in  (12)  contains  an  unlcnov/n 
function  of  t,  and  it  is  convenient  to  combine  this  vilth  other 
functions  of  t  only  and  write  (I2)  in  the  form 
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^  II  dpjd?'  ,  (14) 


v/here 


fit)  =  -  C-  -^-iir^--  +  -  a0[a(t),t]}  . 

(15) 

Integration  of  (l4)  through  the  sphere  at  constant  real  time, 
using  the  boundaiv  conditions  (13)  to  give: 

da 

J.^.V  5/ ‘ 


i3  ,  d  [a(t),t] 


u[a(t),t 


Ja(t)  °P  2'  J.'  ' 

reduces  (l4)  to  an  integral  equation  for  f(t).  Because  ^  is  a 
function  of  r  and  t^  this  space  integral  destroys  the  simplifying 


convolution  property  of  the  integral  equation  over  reduced  time. 


but,  as  shovm  below,  it  can  be  readily  solved  numerically  by 

finite  sum  techniques.  f(t)  having  been  determined,  (12) 

dcJ 

comprises  an  Integral  y  space  integration  from  a  boundary 

gives  <5  ,  and  is  then  determined  from  (4)  . 

I*  y 


Method  of  solution. 

Since  the  functions  which  define  the  problem  —  the 
tempei*ature  field  6(x',t),  the  boundary  radii  a(t)  and  b(t)  and 
the  prescx'lbed  pressures  -fj^(t),  -f^Ct)  --  are  all  defined  in 
terms  of  z'eal  time  t,  and  since  the  benefit  to  be  gained  by  using 
^  as  independent  variable  has  already  been  achieved  through  the 
application  of  convolution  theory  leading  to  (12),  it  is 
advantageous  to  evaluate  the  solution  in  terms  of  the  real  time  t 
Equation  (l^)  then  beccnes 
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dd 

_ r  _ 

SF"  " 


=  -  J  R[C(r,t)-^(r,t»)3|^[f{t‘)  +  r  p3  II  dp]dt».  (17) 
r  o  afi-.t)  °P 


a(t«) 


Integrating  by  parts  to  obtain  a  more  convenient  form  for 
numerical  integration  yields  the  expression 

3jr  =  -  %H(0){f(t)  +  f  p^  ||  dp} 
r  a(t)  ^ 

-  )  {f(fc‘)  +  J  ,  P^  If  dp}  Ip- R{C(r,t)-^(r,t')}dt']. 

o  a(t‘) 

(18) 

As  described  in  [4],  a  suitable  finite  sum  approximation  for  the 
time  Integral  in  (l8)  is  obtained  by  placing  ahead  of  the  Integral 
sign  a  mean  value  of  the  first  bracket  of  the  Integrand.  Equation 
(l8)  can  then  be  written  in  the  form 


where 


Rv  “  R[5(r,t  )  - 

,  ,  3  3e{p,t„) 

h(r,t^)  .  r  p3  5- - dp 

^(tj 


and,  in  order  to  remove  the  constant  6a,  the  functions  have  been 
normalized  with  the  new  given  by  (J^6ci£0^,  where'  E  is-  the 
Instantaneous  Young’s  modulus,  and  0^^  the  majcimum  temperature 
difference  occurring. 
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By  Integrating  the  finite  Evm.  equation  (19)  over  the 
radius,  is  given  in  terms  of  its  value  at  earlier  times; 

hence  the  terms  containing  are  grouped  together,  and  (19) 

is  written  in  the  form 

Rfit  , )  R 

2  =  p  — -TT^  -  f(t  )(D  +  -J)  (22) 

r  r 

v/here 

r^P  =  Ih(r,tj^_^)  +  h(r,t^)  +  ^(^^.1)  ] 


and 


-  (2h(r,t^)  +  f(Vi)Rj,] 
+  Z  [h(r,tj^)  +  h(r,tj^j^)  + 


We  consider  the  case  with  zero  surface  pressure  so  that 
f^(t)  B  fgCt)  =  0,  and  constant  external  radius  b(t)  =  b. 
Integration  of  (22)  over  r  from  a(t)  to  b  then  yields 


pb 

0  =  Pdr  -I-  f(t^_i)G 

i(t) 


Ddr  +  G] 


(23) 


(24) 
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with  the  solution  Imjr.edlately  after  application  of  the  surface 
temperature,  at  t=0  ,  defining  f(o)  =  1.  Substitiatlon  of 
into  (22)  and  Integration  from  the  fixed  outside  boundary  r=b 
determines  according  to 

2dr(r,tn)  =  ^  [f(t^)  -  f  -  JPdr  f(t^)  JDdr. 

(25) 

The  integrals  over  r  v/ere  evaluated  by  using  Simpson's  rule  over 
most  of  the  range,  but  the  trapezoidal  rule  for  a  single  step 
only  and  the  Newton-Cotes  three  interval  expression  adjacent  to 
the  boxindary  r=b  for  odd  numbers  of  intervals  greater  than  one. 
Next  to  the  moving  boundary,  r=a{t),  it  was  necessary  to  use 
unequal  steps  in  r,  and  the  finite  sxim  integration  formula  was 
chosen  to  give  the  area  under  a  parabola  through  the  last  three 
points.  The  choice  of  interval  size  for  both  n  and  t  is  discussed 
in  the  next  section  where  the  solution  is  presented.  The  calcula¬ 
tions  were  carried  out  on  the  IBM-7070  machine  at  the  Brovrn 
University  Computing  Laboratory. 

The  temperature  variation  0(r, t)  was  obtained  as  the 
solution  of  the  heat  conduction  problem,  neglecting  the  coupling 
betv/een  heat  energy  and  mechanical  energy  dissipated  in  visco¬ 
elastic  flow.  For  the  thermal  stress  problem  mentioned  in  the 
Introduction,  we  need  the  temperature  field  in  a  sphere  initially 
at  uniform  temperature,  defined  as  zero,  and  with  the  constant 
temperature  0^^  at  the  ablating  surface.  To  limit  computational 
effort  on  this  particular  problem,  which  is  evaluated  simply  as 


an  exanrple  of  a  class  of  solutions,  a  temperature  field  given  by 
a  relatively  simple  expression  was  achieved  by  considering  the 
cavity  to  be  in  an  infinite  medium  and  by  making  use  of  the 
analogy  betv;een  plane  and  spherically  symmetric  heat  flow 
mentioned  in  [9].  The  resulting  solution  involves  heat  cor-duction 
across  the  radius  r=b,  which  determines  the  outer  bo\andary  condi¬ 
tion  for  the  temperature  distribution  in  the  spherical  shell - 
We  consider  the  temperature  field  for  r  >  a(0)  =  a,  when  the 
variation 

|-=l  +  ^t^  +  Yt  +  |t^  (26) 

is  prescribed  on  the  surface  r=a,  where  p,  y  ^  arbitrary 

constants.  The  locus  r=a(t)  determined  by  0=0^^  can  be  considered 
as  the  ablating  bovindary,  which  can  be  located  from  the  known 
temperature  distribution  corresponding  to  the  boundary  condition 
(26)  at  r=a,  and  zero  initial  temperature  (see  [9]) 

^  =  ^[erfc  X  +  pt^  lerfc  x  +  Yt(4  i^erfc  x) 

^1  ^ 

+  6t^{8  l^erfc  x)  ]  (2?) 

where  x  =  (r-a)/2(Kt)^,  K  being  the  thermal  diffuslvlty.  The 
second  term  with  is  needed  to  ensure  finite  initial  velocity 
for  the  ablating  boundary,  and  the  constants  were  chosen  to 
achieve  approximately  constant  rate  of  ablation.  Figure  1  shows 
the  resulting  curve,  r=a(t),  for  a  =  2.87I,  p=2  and  y=2000, 
compared  v/ith  that  for  constant  velocity  ablation.  For  the 
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nurnGrlcai  solution  the  positions  of  the  ablating  boundary, 
r=a(tj^),  uere  determined  by  solving  (27)  for  x  by  Nev/ton's 
method,  with  the  left  hand  side  replaced  by  unity.  Equation  (27) 
then  determines  the  temperature  field  for  a(tj^)  <  r  ^  h,  and  the 
resulting  distributions  are  shovm  in  Pig.  2  at  four  time  values. 
As  seen  from  Pig.  1,  the  average  velocity  of  ablation  was 
selected  to  complete  the  process,  a(t)=b,  at  t  =  0.7a  /K,  for 
the  outside  boundary  chosen,  b=3a.  Immediately  after  application 
of  the  internal  heat,  at  t=0  ,  the  cavity  surface  is  at  tempera¬ 
ture  with  the  rest  of  the  material  still  at  the  zero  initial 
temperature.  The  temperature  field  does  not  correspond  to  one 
of  the  standard  boundary  conditions  for  heat  conduction,  but  a 
glance  at  the  distributions  adjacent  to  r=3a  suggests  that  the 

temperature  gradient  is  roughly  proportional  to  the  temperature, 

dO 

so  that  the  solution  approximates  the  radiation  condition  ^  = 
-K0.  From  (27)  the  integrals  (21 )  for  h(r,tj^)  can  readily  be 
evaluated. 

As  in  the  problem  discussed  in  [6],  the  thermal  stress 

field  v;as  evaluated  for  polymethyl -methacrylate  material 

behavior  using  measured  tensile  relcixation  data.  We  use  the 

auxiliary  function  R(^)  evaluated  in  [4]  by  numerical  integration 

of  an  Integral  equation  equivalent  to  (ll).  The  shift  function 

9(T)  used  by  Mukl  and  Sternberg  [5]  was  used,  and  also  the  same 

temperature  range,  0^  =  30*^C  with  the  initial  temperature  80°C. 
o 

Taking  K  =  8  cms  /hv,  and  a  =  4  cms,  the  total  time  of  ablation 
becomes  0.7  a  /K  =1.4  hrs .  As  discussed  in  [6],  formally  the 
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solution  applios  for  a  range  of  a  and  K_,  with  constant  a  /X,  but 
since  viscoelastic  characteristics  for  polynethyl -methacrylate 
have  been  used,  for  practical  Interpretation  the  corresponding 
Vcilue  of  K  must  also  be  substituted* 

The  solution 

Using  the  material  data  and  constants  defining  the 
particular  problem  to  be  solved,  time  and  space  steps  for  the 
finite  sirni  integration  procedures  detailed  in  the  previous  section 
were  determined  by  trial .  Increments  were  halved  until  inappre¬ 
ciable  differences  in  the  solution  resulted.  Tlie  full  line  in 
Fig.  3  shows  the  variation  of  f(t)  defined  in  (l5)>  with 
cJj.[a(t),t]  =  0  for  the  case  considered  of  zero  cavity  pressure. 

O 

It  will  be  noted  that,  apart  from  the  factor  -[aCt)]"^,  f(t)  then 
represents  the  excess  of  the  circumferential  strain  at  the 
cavity  surface  over  the  thermal  e:xpansion  there,  and  is  thus  the 
strain  caused  by  stress .  It  was  computed  according  to  (24)  .  As 
detailed  in  the  previous  section,  f(t)  is  the  function  which 
linics  the  stress  values  at  different  radii,  and  constitutes  a 
major  ingredient  of  the  stress  determination.  It  must  therefore 
be  determined  accurately,  and  the  choice  of  time  and  space  steps 
was  assessed  by  their  influence  on  f(t).  The  results  shewn  in 
Fig.  3  were  calculated  for  Ar/a  =  0.05,  and  for  the  early  stages 
KAt/a^  =  0.0125-  The  availability  of  memory  locations  United 
the  time  rai'.ge  over  which  this  fine  time  mesh  could  be  maintained, 
but  increasing  the  step  for  intermediate  times  caused  no  difficulty 
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since  tecgseraturo  gradients,  and  the  resulting  dependent  variables, 
settle  down  to  more  gradual  variation.  When  ablation  Is  almost 
complete,  and  only  a  thin  shell  of  material  remains,  f(t)  changes 
rapidly  and  accuracy  again  demands  a  reduction  in  the  time  step. 
With  only  a  narrow  shell  left,  smaller  increments  in  the  radius 
would  also  be  needed  to  maintain  accuracy,  and  the  completion  of 
the  f(t)  cui*ve  up  to  Kt/a^  =  0.7  was  achieved  by  an  analytical 
expansion  procedure.  That  the  finite  increments  stated  above 
are  small  enough  for  satisfactory  accuracy  Is  indicated  by  the 
comparisons  of  f(t)  values  computed  for  different  finite  differ¬ 
ence  steps  presented  In  Table  I.  With  a  radius  step  of  Ar/a=0.1 
an  undesirable  rovighness  was  apparent  In  the  f(t)  values  due  to 
the  change  in  the  form  of  the  space  Integrals  in  (21 )  and  (24) 
as  the  nxanber  of  radius  increments  spanning  the  ablating  shell 
changed  from  odd  to  even,  and  as  the  size  of  the  interval 
adjacent  to  the  ablating  bovindary  changed.  This  roughness  was 
not  apparent  for  Ar/a  =  .05 . 

Ihe  individual  points  plotted  in  Fig.  3  correspond  to 
the  solution  for  an  elastic  body  with  the  elastic  constants  for 
Instantaneous  loading  of  polyroethyl -methacrylate .  It  is  obtained 
by  simply  substituting  R(?)  =  R(0),  a  constaint,  into  the  visco¬ 
elastic  analysis. 

Computed  values  of  and  <Jq,  evaluated  from  (25)  and 
(4)  are  plotted  in  Pigs.  4  to  7  lor  both  the  viscoelastic  and 
elastic  cases .  Use  of  the  time  and  space  steps  found  satisfactory 
for  determining  f(t)lead  to  satisfactory  accuracy  over  most  of 


Word -1859^'' -6 


l6 


the  radius  range,  but  to  serious  inaccuracies  adjacent  to  the 
ablating  boundary  in  the  viscoelastic  case.  Those  occur  because 
the  temperature  rise  at  the  ablating  surface  results  in  acceler¬ 
ated  stress  relaxation  due  to  the  sensitivity  of  the  scale 
■factor  (p(T)  to  temperature.  Taking  the  initial  temperature  as 
the  base  temperature,  9  is  lO'^*  at  the  ablating  surface,  so 
that  relaxation  times  are  all  reduced  by  this  factor  with  the 
resulting  rapid  stress  relaxation  shown.  In  order  to  improve 
accuracy  in  the  region  of  the  ablating  bovindary,  integration  of 
(22)  for  was  carried  out  from  the  inner  boundary  a(t)  using 
time  steps  of  .00125  a  /K,  and  space  steps  of  Ar/a  =  .0125. 
Corresponding  values  of  f(t)  were  deduced  by  interpolation. 
Results  of  these  calculations  are  shown  by  the  circled  points  in 
Pig.  4,  and  they  are  seen  to  veer  smoothly;  into  the  values 
calculated  for  the  larger  mesh  away  from  the  inner  boundary. 
Points  corresponding  to  the  latter  are  indicated  by  crosses  and 
are  seen  to  be  appreciably  in  error  near  the  boundary.  Similar 
corrections  were  needed  at  later  values  of  the  time,  though 
somewhat  larger  time  steps  could  be  tolerated,  and  were 
necessitated  by  machine  storage  limitations,  unless  more 
elaborate  programming  techniques  were  to  be  utilised.  As. men¬ 
tioned  above,  smaller  time  and  radius  steps  were  needed  to 
maintain  accuracy  near  the  end  of  the  process  when  only  a  thin 
spherical  shell  remained,  and  these  could  be  incorporated  since 
only  a  small  radius  range  occurred . 
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Study  of  the  resulto  depicted  in  Pigs.  3  to  7  reveal 
interesting  features  of  the  solutions.  The  stress  distributions 
show  that  the  Influence  of  viscoelasticity  is  large  but  localized 
in  the  heated  region  near  the  ablating  boundary,  in  confonnity 
with  the  sharp  Influence  of  temperature  on  viscoelasticity 
through  the  scale  factor  9(T)  .  The  main  body  of  the  material 
thus  responds  essentially  elastically.  The  thermal  expansion 
near  the  ablating  surface  is  resisted  by  the  constraint  of  these 
outer  layers,  and  this  leads  to  the  high  peak  of  circumferential 
compression.  Immediately  on  application  of  the  cavity  surface 
temperature  0^^,  at  t=o'  ,  the  strain  is  completely  constrained  to 
remain  zero,  and  the  cavity  surface  la  subjected  to  a  surface 
compression  of  magnitude  0.2564,  as  discussed  in  [6]  for  the 
solid  sphere.  Qualitatively,  constraint  of  this  nature  remains 
until  ablation  is  almost  complete  as  indicated  by  the  variation 
of  f(t),  Vfhlch,  as  mentioned  above.  Is  a  function  of  the  cir¬ 
cumferential  expansion  at  the  cavity  surface  due  to  stress  alone . 
For  quantitative  assessment,  the  factor  -[a(t)]^  must  be  taken 
into  account.  The  rapid  decrease  in  f(t)  ‘shov/n  in  Fig.  3  when 
ablation  is  nearly  complete  represents  freedom  for  thermal 
expansion  \/hen  the  temperature  rise  spreads  to  the  outer  boundary 
and  a  consequent  decrease  in  magnitude  of  the  negative  strain  due 
to  the  compressive  stress.  In  the  elastic  case  the  strain  just 
prior  to  complete  ablation  is  entirely  due  to  thermal  expansion, 
since  stresses  reduce  to  zero  with  no  surface  pressure,  acting, 
and  hence  f(tjj)  »  0,  where  a(t^)  =  b.  For  the  viscoelastic  case. 
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the  outer  ring  has  been  stressed  In  tension  throughout  the 
process,  and  so  a  positive  residual  strain  due  to  viscoelastic 
flov;  is  superposed  on  the  thermal  strain,  and  f(tjj)  becomes 
negative  as  shown. 

An  Indication  of  the  effectiveness  of  the  constraint 
of  the  outer  part  of  the  sphere  is  provided  by  comparison  with 
the  elastic  solution  for  a  cavity  in  an  infinite  solid.  In  this 
case,  irrespective  of  the  form  of  the  spherically  symmetrical 
temperature  distribution,  Sternberg  [lO]  has  shovm  that  the 
cavity  surface  displacement  remains  zero.  Because  a  quasi-static 
elastic  solution  depends  only  on  the  current  geometry  and 
tractions,  it  follows  that  u[a(t),t]  =  0  for  the  elastic  solution 
corresponding  to  our  ablation  problem  but  with  the  outer  radius 
b  increased  to  Infinity.  Substituting  this  into  the  expression 
for  f(t)  yields  the  broken  curve  shov/n  in  Fig.  3*  It  lies  close 
to  the  finite  sphere  solution  at  small  times,  aind  the  latter 
gradually  falls  below  it  due  to  the  positive  radial  displacement 
permitted  by  the  removal  of  restraint  forces  at  the  outer  bound¬ 
ary,  It  is  seen  that  this  release  is  quite  small  until  appreci¬ 
able  ablation  has  taken  place.  The  reduced  peak  of  circumferen¬ 
tial  compressive  stress  in  the  viscoelastic  case  reduces  the 
tendency  for  the  outer  shell  to  expand,  hence  reduces  the 
internal  displacement,  and  so  the  viscoelastic  cui^ve  in  Fig.  3 
lies  between  the  elastic  values  and  that  for  the  infinite  elastic 


sphere . 
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The  stress  distributions  for  the  ablating  hollow 
sphere  presented  here  contrast  v;lth  the  solid  sphere  case 
described  in  [6].  In  the  present  case  the  thermal  stress  field 
dies  out  much  less  rapidly  than  for  the  solid  sphere.  For  the 
solid  sphere  the  heat  from  the  ablating  outer  surface  heats  up 
the  entire  sphere  and  the  temperature  approaches  the  uniform 
value  With  decreasing  temperature  gradients,  the  thermal 

stresses,  even  in  the  elastic  case,  decrease  rapidly  with  time. 
The  condition  of  heat  transfer  across  the  outer  boundary  in  the 
present  case  causes  retention  of  high  temperature  gradients  as 
is  clear  from  Pig.  2,  and  this  leads  to  persistence  of  thermal 
stresses  in  both  the  elastic  and  viscoelastic  cases.  Thermal 
insulation  of  the  outer  boundary  would  reduce  these. 

As  in  [6]  the  dimensions  of  the  sphere  (a=4  cms, 
b=12  cms)  and  the  temperature  field  (total  ablation  time  1.4  hrs.) 
were  chosen  to  illustrate  a  significant  influence  of  thermo- 
viscoelasticity.  A  shorter  ablation  time  would  reduce  the  influ¬ 
ence  of  viscoelasticity  foi'  the  sa.me  temperature  range  and  scale 
factor  9(T).  The  modest  temperature  range  80®C  to  110°C  was 
chosen  for  comparison  v/ith  previous  solutions  ([5]  and  [6]),  and 
the  effect  of  temperature  on  viscoelasticity  is  so  great,  that 
an  Increase  in  the  upper  temperature  of  the  order  100°C  could 
increase  9  by  several  factors  of  10,  and  yield  a  marked  visco¬ 
elastic  effect  even  for  ablation  times  of  the  order  of  seconds 
In  the  same  material .  Thus  quantitative  estimates  of  dominant 
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relaxation  times  and  the  variation  of  the  scale  factor  cp(T)  are 
needed  to  assess  the  expected  Influence  of  viscoelasticity  on  a 
thermal  stress  problem. 

This  solution  has  been  presented  to  illustrate  the 
efficacy  of  the  Integral  equation  approach  to  a  thermo-visco¬ 
elastic  stress -analysis  problem.  Humerlcal  solution  can  be 
readily  processed  on  an  electronic  digital  computer,  and  this 
enables  experimentally  measured  viscoelastic  characteristic  to 
be  directly  Introduced  into  the  analysis.  This  solution  goes 
beyond  those  previously  published  in  this  area,  since  the 
analysis  reduces  to  a  integral  equation  spanning  both  the  space 
and  time  variables . 


&B9 

.0125 

.025 

— 

.05 

.05 

1 .592 

f(.i) 

2.450 

2.442 

2.431 

wasm 

-- 

r{A) 

-- 

— 

8.350 

.10 

f(.05) 

1.559 

1.562 

1.593 

f(.i) 

2.415 

2.412 

2.411 

f(.2) 

4.580 

4.563 

4.526 

f(.4) 

— 

8.369 

8.321 

Table  I.  f(t)  computed  using  different  time  and  space  steps. 
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FIG.  2  TEMPERATURE  PROFILES  AT  VARIOUS  TIMES 
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FIG,  3  VARIATION  OF  f(t)  =  -  |}  WITH  ^ 
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FI6.4  STRESS  DISTRIBUTION  FOR  SPHERE  b/a  =  3  AT  TIME  .05  oVk 
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FIG.  5  STRESS  DISTRIBUTION  FOR  SPHERE  bAi=3  AT  TIME  .20  aVK 
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FIG.  6  STRESS  DISTRIBUTION  FOR  SPHERE  b/o  =  3  AT  TIME  .40  o'/K 
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FIG.7  STRESS  DISTRIBUTION  FOR  SPHERE  b/a=3  AT  TIME  .60  a*/K 


